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G\ ; Abstract 

On 

We examine whether the self interaction correction method by Harrison, 
►_> ■ which does not introduce the spherical single particle density approximation 

\Q . to energy functional, can be applied to Na clusters. We show that it does 

not work well, especially, for large clusters, though it works well for atomic 
systems. We suggest that it is better to apply this method only to the Hartree 
term. We also show that the effects of non-diagonal Lagrange multiplier 
-y-, i originating from the orthonormality of single particle orbitals are negligible. 
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I. INTRODUCTION 

The local density approximation (LDA) provides a powerful practical technique to apply 
the Kohn-Sham framework |IJ to interacting many body problems. A problem of this method 
is unphysical self interaction, i.e. the interaction of a particle with itself. Perdew and Zunger 
0H proposed a prescription to remedy this shortcoming, which has been used for atoms 
0-P, molecules @J^|, bulk systems |||| and also metal clusters [p^H l2|]. Though there still 



remain some effects of self interaction, the major part of the problem is removed in this 
method. Compared with the Hartree-Fock theory, which is free from the self interaction 
problem, the local density approximation with the self interaction correction (SIC) has 
advantages such as, i) exchange and correlation energies can be relatively easily handled 
in the same manner, ii) the resultant single particle energies well approximate the physical 
removal energy from each orbit, iii) the numerical calculation is much lighter, especially for 
three dimensional calculations. 

A characteristic feature of the SIC method of Perdew and Zunger is that the energy 
functional depends not only on the total density, but explicitly also on the density of each 
single particle orbital. In almost all calculations for closed shell atoms and metal clusters, the 
single particle densities in the energy functional are substituted by the spherically averaged 
densities. The central single particle potential is then deduced by taking functional derivative 
of the resultant energy functional with respect to the spherically averaged single particle 
density. Following more closely the original idea of Perdew and Zunger, on the other hand, 
Harrison j|,|5| proposed a method of using the original single particle densities without 
introducing spherical averaging. Since the energy functional is not invariant under unitary 
transformation of single particle orbitals, Harrison represented the single particle orbitals 
by either spherical harmonics or Cartesian basis as two choices. 

Though the method by Harrison works well for atoms 0||, it has not been tested for 
metal clusters. We address this question in this paper by taking Na clusters as an example. 
We show that it does not work well, especially for large clusters which have single particle 
orbitals with large angular momentum. We confine our study to the exchange energy without 
referring to the correlation energy in order to make the argument clear and compare the 
results with those of Hartree-Fock calculations. We show that a better agreement with the 
Hartree-Fock calculations is obtained if one applies Harrison's method only to the Hartree 
term. 

In addition to the validity of Harrison's method, we discuss in this paper the problem 
of non-diagonal Lagrange multipliers originating from the orthonormality of single particle 
orbitals in the self interaction correction method of Perdew and Zunger. We show that the 
non-diagonal property of the Lagrange multipliers introduces a negligible effect to single 
particle energies as well as the total energy. 

The paper is organized as follows. In Sec.|TI| the SIC method of Perdew and Zunger and 



Harrison's approach are briefly explained. In Sec . [III A] the exchange energies calculated by 
several SIC methods are compared, and the effect of non-diagonal Lagrange multipliers is 
discussed in Sec.[TirB]. Summary and conclusion are given in Sec .[TV]. 



II. SIC FORMALISM AND HARRISON'S METHOD 

The total energy in the self interaction corrected local density approximation (SIC-LDA) 
by Perdew and Zunger 0|| is expressed as 

-■sic 



where 



E TOT = T + E ext + E H + E x , (1) 

T = -^E/^^WV 2 ^(r), (2) 
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E cxt [p] = J d 3 r v cxt (r)p(r) , (3) 

E M = y*r J &!&*£, (4) 
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ET = E^[ Ph Pi ] - £ {E H [ Pi ] + E^[ Pu 0]} , (5) 
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E^ A \p hPi ] = -\ (£) 3 E /^ ^W 4/3 • ( 6 ) 



Here all quantities are in Hartree atomic units, i.e. m = e 2 = h = 1. As mentioned in 
the introduction, the correlation energy has been neglected. We take the effects of ions into 
account in the spherical jellium model. The external potential is then given by 

Vext(r) = { -Z/(2R icll ) {3 - (r/i? Jcll ) 2 } r < i? jell 
1 -Z/r r > i? je ii , 

where the jellium radius i?j e ii is related to the number of atoms in the cluster Z by i?j e u = 
r s Z 1//3 , r s being the bulk Wigner-Seitz radius which is 4 a.u. for Na. 
The Euler equation under the orthonormality condition 



J E tot + J2 e* (tij ~ f d 3 r ^*(r)^(r; 



results in the following coupled equations for the single particle wave functions 
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The Lagrange multiplier e^ becomes non-diagonal because of the orbital dependence of the 
self interaction corrected exchange potential. In the following, we approximate it by the 
diagonal components and solve the following non-coupled equations [f| 



-^V 2 + ^xt(r) + /rfV. 
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We discuss the validity of this approximation in Sec. |.l nB| for Na clusters. 

As we see in Eq.(|5]) the SIC method of Perdew and Zunger is characteristic in that the 
total energy functional depends explicitly on the individual orbital density. A consequence 
is that it loses invariance under the unitary transformation of single particle orbitals. 

Another problem is that the numerical load is heavy, because one has to solve three- 
dimensional equations instead of the one- dimensional equations for the radial motion of 
electrons even for closed shell atoms and metal clusters. In applying this formalism to those 
systems, one usually replaces p%{r) in the curly brackets in Eqs.(|[) and ([U]) by the spherically 
averaged orbital density given by Pi(r) 



Pi[r) 



^Jdr Pi (r). 
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This prescription certainly reduces the numerical load, because the resultant single particle 
potentials become central potentials. However, it does not optimize the SIC following the 
original scheme of Perdew and Zunger. 

Harrison proposed an alternative procedure, which avoids replacing the single particle 
densities in the energy functional by spherically averaged ones. As mentioned in the intro- 
duction, he expressed them in the spherical harmonic basis 
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or in Cartesian basis 
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He noticed that the resultant energy functional can be expressed in terms of the spherically 
averaged orbital density after the integration over angle, 



ET = E^[ Ph Pi ] - £ 2(2/ + 1) (f; E%\p. 

n,l Kk=0 

E H l [pnl] = 2° k H J dr J dr ' KK r )| 2 \ U nl{r')\ 



I 7T.LDAr~ 



c l x E x 



[Pnl,0] 



2 < 



(15) 
(16) 



k,l 



The coefficients c# and c l x for each of the spherical harmonic and Cartesian representations 
are listed in Tables |I| and ||. Harrison then calculated the corresponding central potential 
for each set of quantum nembers n, I by taking the functional derivative with respect to the 
spherically averaged orbital density p n i(r) 



spherical Cartesian 

harmonics 



1 2 0.0800 0.1600 

2 2 0.0571 0.0816 

2 4 0.0317 0.0816 

3 2 0.0533 0.0686 
3 4 0.0202 0.0346 

3 6 0.0179 0.0538 

4 2 0.0519 0.0632 
4 4 0.0180 0.0269 
4 6 0.0108 0.0208 
4 8 0.0119 0.0398 

TABLE I. The coefficients c^ for the Hartree term. 



spherical Cartesian 

harmonics 



1 1.0937 1.1800 

2 1.1293 1.2384 

3 1.1508 1.2708 

4 1.1659 1.2925 

TABLE II. The coefficients c x for the exchange term. 
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SIC(nl), , &E X [p rui 

v x (r) = A ~ , v • (17) 

Though this method restricts the variational space smaller than that in the original 
scheme of Perdew and Zunger since it presumes a spherically symmetric potential from the 
beginning, Harrison showed that his method still improves both the exchange energy and 
the total energy for atoms compared with the simple procedure where Pi(r) is replaced by 
the spherically averaged orbital density. 

It is an interesting question to see whether Harrison's method can be applied to metal 
clusters. A simple minded consideration would suggest that Harrison's method becomes 
more powerful in metal clusters. This is because Harrison's treatment should have a large 
effect on high angular momentum orbitals which play more important roles in metal clusters 
than in atoms where the main contribution to the SIC originates from the Is state 0. In 
the next section, we apply Harrison's method to Na clusters, and show that it does not work 
well contrary to the simple expectation. 



Nag Na2o Na4o Nag2 

HF -2.95 -2.95 -2.95 -3.03 

LDAX -2.68 -2.78 -2.84 -2.93 

SA-SICX -2.93 -2.91 -2.91 -2.95 

SH-SICX -2.87 -2.84 -2.83 -2.87 

(H only) -3.00 -2.96 -2.96 -3.00 

C-SICX -2.82 -2.78 -2.77 -2.80 

(H only) -3.07 -3.01 -3.00 -3.02 

TABLE III. Exchange energy per electron in eV. SA-SICX, SH-SICX and C-SICX are the 
abbreviations of the SICX calculated with spherically averaged, spherical harmonic and Cartesian 
orbital densities. H only means that only the Hartree term has been calculated by using either 
spherical harmonic or Cartesian bases (see text). 



III. RESULTS AND DISCUSSION 

A. Exchange and total energy 

We compare in Table [III] the exchange energy per electron for Na clusters calculated 
by several methods. In this table and in what follows, the abbreviations SA-SICX, SH- 
SICX and C-SICX stand for the SICX calculation using the spherically averaged, spherical 
harmonic and Cartesian orbital densities, respectively. The difference between the HF and 
the other calculations represents the error of each method since the HF calculation provides 
the exact exchange energy. Strictly speaking, the exact exchange energy in the Kohn-Sham 
formalism is the one given by the optimized effective potential method |]13||. It is known, 



however, that it is nearly the same as that given by the HF calculation for atoms [0,0 . 

We first compare the results of HF, LDAX and SA-SICX for four different Na clusters. 
The order of the estimated exchange energies, which are negative, is LDAX > SA-SICX > 
HF irrespectively of the size of the Na cluster. This is different from the order for atoms, 
where LDAX > HF > SA-SICX [|. The deviation of the result of LDAX from that of HF 
gets smaller with the size of the Na cluster, while that of SA-SICX is almost constant. 

We next compare the results of HF and three SIC methods. The absolute value of the 
exchange energy calculated by using the spherical harmonic and Cartesian orbital densities 
becomes smaller than that estimated by the SA-SICX. Including the negative sign, the order 
is C-SICX > SH-SICX > SA-SICX. Consequently, the deviation of the results of SH-SICX 
and C-SICX from that of the HF gets larger than that for SA-SICX. These two calculations 
are even worse than LDAX for large systems. Their deviation from the HF calculation gets 
larger with increasing size of the cluster. 

For atoms, the order of the exchange energies C-SICX > SH-SICX > SA-SICX is the 
same as that for Na clusters. However, the result of C-SICX is still below that of HF. This 
means that both of SH-SICX and C-SICX have smaller deviation from the HF result than 
SA-SICX and provide better prescriptions for atoms. 
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FIG. 1. Each orbital contribution to the self interaction energy for Na4o. SA, SH and C 
mean the SICX with spherically averaged, spherical harmonic and Cartesian orbital densities, 
respectively. 



We now investigate the reason why the SH-SICX and C-SICX methods are inferior to 
SA-SICX for Na clusters. Since the self-consistent density in the three SICX methods is 
almost the same, the first term of Eq.flTSp is nearly the same for the three methods. The 
major difference among them is therefore associated with the second term of Eq.([TJ|), which 
consists of contributions from orbitals with various n, I quantum numbers. In order to see 
the physics clearly, we consider the following self interaction averaged over the azimuthal 
quantum number for each set of quantum numbers n, I, 



21 + 



— 7 /I \EH\pinm] + E x [pnlm, 0] j . 



(18) 



Its value evaluated by inserting the orbital densities given by Eqs.(|l2|), (|13|) and (|I4D is 
compared in Fig.|T| for the occupied orbitals in Na^. 

The three methods agree well for the s orbitals. The small differences originate from 
a subtle difference of the total densities in the three methods. On the other hand, a large 
difference appears among them for the finite angular momentum states, i.e. for the p, d and 
/ orbitals. Two specific features can be remarked, i) always SA-SICX > SH-SICX > C- 
SICX, and ii) the higher the orbital angular momentum is, the larger the differences among 
them are. 

The key point to understand these features is the degree of localization of single particle 
density adopted in the three SICX methods. Clearly they differ in their angular treatment 
of the orbital density. C-SICX has the most distinct localization in angle, while SA-SICX 
has, of course, no angular localization. Symborically, we express this situation of the degree 
of localization as C-SICX > SH-SICX > SA-SICX. This difference gets more prominent for 
higher angular momentum. 

The problem is that the angular localization property affects the self Hartree term, i. e. 
the first term in the curly brackets in Eq.([5J), and the self exchange terms, i.e. the second 
term, in different way. Since the self Hartree term uses the exact long range Coulomb 
interaction between electrons, it does not depend so much on the density profile of electrons 
including the angular localization. On the other hand, the self exchange term is evaluated 
based on the LDA. This is equivalent to assuming a short range S interaction, which leads to 





Na 8 


Na 20 


Na 40 


Na 92 


HF 


-140.9 


-626.2 


-1957.6 


-7766.7 


LDAX 


-138.9(+1.42%) 


-623.0(+0.51%) 


-1953.3(+0.22%) 


-7758.1(40.11%) 


SA-SICX 


-140.8(+0.07%) 


-625.7(+0.08%) 


-1956.4(+0.06%) 


-7761.1(40.07%) 


SH-SICX 


-140.4(+0.35%) 


-624.4(+0.29%) 


-1953.5(+0.21%) 


-7754.0(40.16%) 


(H only) 


-141.3(-0.28%) 


-626.8(-0.10%) 


-1958.4(-0.04%) 


-7765.1(40.02%) 


C-SICX 


-140.1(+0.57%) 


-623.2(+0.48%) 


-1951.0(4-0.34%) 


-7747.7(40.24%) 


(H only) 


-141.8(-0.64%) 


-627.7(-0.24%) 


-1960.0(-0.12%) 


-7767.5(-0.01%) 



TABLE IV. Total energy for Na clusters in eV. The relative error of each method (the difference 
between the values in each method and HF divided by the HF value) is shown in parentheses. 



a strong sensitivity of the exchange energy on the details of the density profile of electrons. 
The more the orbital localizes, the larger the absolute value of the self exchange energy is. 
In C-SICX, which has the most prominent angular localization, the negative self exchange 
energy increases with angular momentum and eventually overwhelms the self Hartree term 
resulting in the sign change of the total self interaction energy. This can be clearly seen 
in Fig.|l]. A less prominent, but a similar, trend can be seen also for SH-SICX. However, 
this strong sensitivity of the self correction energy to the angular localization of electronic 
density may be unphysical because the LDA for the exchange term cannot be justified for 
the localized density in SH-SICX and C-SICX. The use of spherically averaged density in 
the SA-SICX moderates to some extent the overamplification of the angular localization 
dependendence leading to a better approximation. 

The localization of radial wave functions in atoms is very different from that in Na 
clusters when one compares the states with the same nodal quantum number. For example, 
Is and lp orbitals have very similar radial distibution in Na clusters, while the latter (2p in 
the atomic notation) is much more extended than the former in atoms. Consequently, the 
self interaction correction mainly originates from the most localized Is state || in atoms. 
Moreover, only small angular momentum orbitals appear in atoms. Thus, the problem stated 
above, i.e. the unphysical strong angular localization dependence does not cause so serious 
trouble in atoms. 

To remedy the problem for Na clusters, we propose to utilize spherical harmonic and 
Cartesian orbital densities only for the self Hartree term. This is a natural consequence of 
the considerations mentioned above. The exchange energy calculated in this way is given in 
the lower side in the rows for SH-SICX and C-SICX in Table |II. They are designated by 
the label "(H only)". All of the H-only exchange energies better agree with that in HF than 
the corresponding results of the SH-SICX and C-SICX which have been calculated using 
the spherical harmonic or Cartesian orbital densities both for the self Hartree and exchange 
terms (the upper side in each row). Especially, H-only SH-SICX is superior to the SA-SICX 
for most systems. We remark that a similar improvement has been obtained by the H-only 
method concerning the total energy. 



B. Effects of non-diagonal Lagrange multipliers 

Before we conclude the paper we comment on the effects of off-diagonal Lagrange mul- 
tipliers in Eq.(|]). We performed calculations by keeping off-diagonal Lagrange multipliers 
for Na 20 and Na 40 with SA-SICX, SH-SICX and C-SICX. We found that the differences of 
the total energy and single particle energy obtained by self-consistently solving Eq.(|9D and 
Eq.fllTl) are less than 10 _1 eV and 10~ 2 eV, respectively. We therefore conjecture that one 
can safely ignore the off-diagonal components of the Lagrange multipliers for Na clusters. 
This result is consistent with that in Ref. H for atoms. 



IV. SUMMARY AND CONCLUSION 

We have calculated the self interaction corrected exchange energy for Na clusters by using 
spherically averaged, spherical harmonic and Cartesian orbital densities. We found that 
both calculations using the spherical harmonic and Cartesian orbital densities deviate from 
the HF results more than the calculations with spherically averaged orbital densities. The 
deviation is especially large for systems with large angular momentum orbitals. We attribute 
this problem to the LDA to evaluate the self exchange term. From this consideration, we 
propose to use the spherical harmonic and Cartesian orbital deinsities only for the self 
Hartree term and to use spherically averaged orbital densities for the self exchange term. 
We have shown that this treatment improves indeed the exchange energy to well reproduce 
the results of HF calculations. We expect that the remaining errors can be diminished by 
the SIC using GGA (generalized gradient approximation) [|l5H i8| . 
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